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Abstract 

We study the spectrum of the integrable open XXX Heisenberg spin chain sub- 
ject to non-diagonal boundary magnetic fields. The spectral problem for this 
model can be formulated in terms of functional equations obtained by separation 
of variables or, equivalently, from the fusion of transfer matrices. For generic 
boundary conditions the eigenvalues cannot be obtained from the solution of 
finitely many algebraic Bethe equations. Based on careful finite size studies of 
the analytic properties of the underlying hierarchy of transfer matrices we devise 
two approaches to analyze the functional equations. First we introduce a trun- 
cation method leading to Bethe type equations determining the energy spectrum 
of the spin chain. In a second approach the hierarchy of functional equations is 
mapped to an infinite system of non-linear integral equations of TBA type. The 
two schemes have complementary ranges of applicability and facilitate an efficient 
numerical analysis for a wide range of boundary parameters. Some data are pre- 
sented on the finite size corrections to the energy of the state which evolves into 
the antiferromagnetic ground state in the limit of parallel boundary fields. 

PACS: 02. 30. Ik, 75.10.Pq 
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1 Introduction 

The solution of the spectral problem for integrable models constructed within the framework 
of the Quantum Inverse Scattering method (QISM) is facilitated by a collection of various 
very powerful tools, commonly subsumed as Bethe ansatz methods. In Bethe's original work 
on the spin-1/2 Heisenberg chain pE] eigenf unctions were found by means of an ansatz for 
solutions of Schrodinger's equation in terms of scattering states of m magnons created from 
the completely polarized state as a reference which led to a system of coupled algebraic 
equations for m parameters. The roots of these Bethe equations correspond to particular 
eigenstates of the model in the m-magnon sector. Within the QISM this approach was put on 
an algebraic basis. For representations of a Yang-Baxter algebra the hamiltonian is identified 
as one member of a family of commuting operators generated by the transfer matrix. The 
eigenstates of this transfer matrix are obtained by the action of creation operators on the 
ferromagnetic reference state. The creation operators are functions of the roots of the Bethe 
equations. 

This need for an eigenstate which is sufficiently simple to guess, however, limits the use 
of both the coordinate and the algebraic Bethe ansatz. For other cases, e.g. for models 
where the underlying symmetry is realized in a non-compact way without a highest weight 
state [31 [23] or situations where the U(l) symmetry of the bulk is broken by interactions 
or boundary terms [18j[T9] different approaches are needed: within the class of functional 
Bethe ansatz methods, e.g. Baxter's method of commuting transfer matrices, Separation of 
Variables (SoV), or the fusion algebra [S1II11ES], relations between elements of the Yang- 
Baxter algebra are derived which eventually lead to functional equations for the eigenvalues 
of the transfer matrix. These relations have to be solved based on the analytic properties 
of the eigenvalues, e.g. distribution of zeros and poles, and their asymptotic behaviour. 
In certain cases this approach may lead to Bethe equations similar to those found within 
the coordinate or algebraic Bethe ansatz. In some cases where this is not possible, the 
functional relations have been brought to a form which allows to express them in terms of 
non-linear integral equations. From the analysis of these equations it has been possible to 
gain important insights into the properties of the ground state and low lying excitations as 
well as the thermodynamics for certain models, see e.g. |12j . 

In this paper we employ some of these functional methods to study the spectrum of the 
isotropic spin-1/2 Heisenberg chain with open boundary conditions where the symmetry of 
the bulk is broken due to non-parallel boundary magnetic fields. The model is given by the 
hamiltonian 
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where a", a = x, y, z denote the Pauli matrices acting on the space of states of a spin-1/2 at 
site j and a^, /3 ± , 9^ parametrize the boundary fields acting on sites L and 1 respectively. 
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The hamiltonian (11. ip is obtained from a transfer matrix based on a representation of 
Sklyanin's reflection algebra [24] which extends the QISM to systems with open boundary 
conditions. Thereby the integrability of the model is established. At the same time, however, 
the fact that the ferromagnetically polarized state with all spins up is not an eigenstate of 
%xxx for generic boundary fields prevents the application of the coordinate or algebraic Bethe 
ansatz to the solution of the spectral problem. In a previous paper [9] we have used Sklyanin's 
Separation of Variables method to address this problem and have derived difference equations, 
so-called TQ-eqyidXioTis, which are satisfied by all eigenvalues of the transfer matrix. To access 
properties of the system in the thermodynamic limit L — Y oo, however, the characterization 
of the Q-functions appearing in these equations was incomplete. This is a familiar situation 
arising in the actual computation of spectral properties: to extract useful information from 
the Bethe equations one always needs some additional insights into the behaviour of their 
solutions as the system sizes varies, in particular for those corresponding to states with low 
energies or to the equilibrium state at finite temperature. This requirement appears to limit 
several other attempts to study spin chains with non-diagonal boundary terms based on e.g. 
representations of a g-Onsager algebra [2] or on an alternative use of the Yang-Baxter algebra 
due to Galleas [ID] , Here we tackle this difficulty by identifying the analytical properties of 
the objects appearing in the TQ-equations and the related fusion hierarchy from finite size 
studies. 

Our paper is organized as follows: in Section [2] we present a brief account of the construc- 
tion of integrable boundary conditions within the QISM. We recall previous results obtained 
from SoV and by fusion of transfer matrices with higher dimensional auxiliary spaces. The 
spectral parameter in the TQ-equations derived within the SoV approach is restricted to a 
finite lattice of points which prevents the use of the functional methods for their solution. 
On the other hand, within the fusion procedure an equivalent TQ-equation with continuous 
arguments is obtained assuming that a certain limit of the fused transfer matrices exists for 
infinite dimensional auxiliary space (see also [21j[29]). Based on this equivalence we perform 
a finite size study of the analytical properties of the transfer matrix eigenvalues and of the 
Q-functions in Section [3J This allows one to derive a finite system of 'truncated' Bethe equa- 
tions which we solve numerically for selected eigenstates of (jl.ip . In Section 0] we use our 
finite size data to rewrite the fusion hierarchy in terms of non-linear integral equations of 
TBA type whose solution determines one selected eigenstate of (II. ip . The paper ends with 
a summary of our results and some concluding remarks. 

2 Integrable Boundary Conditions 

The construction of integrable systems involving boundaries within the QISM was initiated 
by Sklyanin |24j. It is valid for a general class of integrable systems characterized by an R- 
matrix of difference form R(X, /i) = R(\ — /i)G End(V" <g) V) which satisfies the Yang-Baxter 
equation 

R 12 (\- f j,)R 13 (\-v)R 23 (n-v) =R 23 (v-v)R 13 (\-v)R 12 {\-fj,) . (2.1) 

The indices of Rjk denote the embedding where R acts non-trivially on the tensor product of 
vector spaces V\ <g> V 2 <8> V3. For the hamiltonian (11.11) we need the well-known 6- vertex model 
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solution 

/ a ( A ) \ a(A) = A + i 

b(X) c(A) 

c(A) b(X) 

\ a(X)J 



R(X) 



6(A) = A (2.2) 
c(A) = i 



of the Yang-Baxter equation ( 12.11) with Vj = C 2 . Each solution -R(A) fixes the structure 
constants of the related Yang-Baxter algebra 

R 12 (X - /i)Ti(A)T 2 (//) = T 2 ( f i)T 1 (X)R 12 (X - /z) (2.3) 

with generators T%X), a,f3 = 1,2. T X (A) = T(A) ®J, T 2 (A) = Ig)T(A) are the embeddings 
of the monodromy matrix T(A) in the product of auxiliary spaces V± <8> V2, they satisfy the 
inversion formula 

r w=($o «w) • rl{x ^w)hwr Tt{x -' l)aV ■ {2A) 

The scalar factor (d q T)(X) = A(X + i/2)D(X - i/2) - B(X + i/2)C(A - i/2) is the central 
element of the Yang-Baxter algebra and is known as quantum determinant. The superscript 
t denotes a transposition and can be extended to the j'th auxiliary space by tj . 

Sklyanin's construction of open spin chains is based on the representations of two algebras 
IA + {X) and IA~(X) defined by the relations 

R 12 {X - fi)U^{X)R 12 {X + fi- i)W 2 ~(/i) = U 2 {p)R 12 {X + // - i)Wf {X)R 12 (X - //) (2.5) 

R 12 { fJl - X)U+ l1 (A)i? 12 (-A -fi- i)U 2 +t2 (/i) = U 2 +t2 (fj,)R 12 {-X -fi- i)U+ fl (A)i? 12 (/i - A) (2.6) 

We shall call U + (X) and U~(X) right and left reflection algebras respectively. Considering 
their product in auxiliary space the trace 

r(A) =trW + (A)W~(A) (2.7) 

defines the transfer matrix as the central object under consideration because it generates with 
[t(A), r(fi)] = a commuting family of operators which can be simultaneously diagonalized. 

The explicit construction of integrable open boundary conditions for models arising from 
the Yang-Baxter algebra with the i?-matrix (12.21) starts with the 2x2 matrix 



K( \ i \ _ ± f a± ch + A sh f}^ Xe e , , 

l ' ;_ a±ch/3±l Ae- e± a± ch (5 ± - X sh f3 ± ' 



independently found by [7] ITT], considered in the rational limit [15] and used here in a 
parametrization introduced by Nepomechie [20]. It constitutes the known c- number rep- 
resentations K + (X) = K(X + i/2, +) and K~(X) = K(X — i/2, — ) of the reflection algebras 
( 12. 5p . ( 12 .6p with the properties ti K(X, ±) = 2 and K(0, ±) = I 2 . Following Sklyanin we 
choose 

U-{X) = T(A - i/2)K(X - i/2, -)o<T*(-A - i/2)a* = ^~J^j (2.9) 
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as representation incorporating the inversion formula (12. 4p . This leads to an explicit represen- 
tation of the transfer matrix with normalization condition r(i/2) = (d q T)(— i/2) responsible 
for an additional factor of 1/2 in a similar decomposition of the right reflection algebra 

«+(A) = \l<(X + i/2, +) = . (2.10) 

Then the hamiltonian (II. ip is connected to the transfer matrix by Hxxx = i<91n r(i/2). 



2.1 Quantum Determinants 

Analogously to the quantum determinant of the Yang-Baxter algebra there exists similar 
objects for the reflection algebras. Quantum determinants play a crucial role when applying 
functional methods to solve the spectral problem. For the left reflection algebra it is defined 
according to |24j reading with the projector P{~ 2 onto the singlet in V\ <8> V2 

(A"W)(A) = tr 12 Pgj{ (A - i/2)i? 12 (2A - i)W 2 "(A + i/2) . (2.11) 

To express (A~W)(A) in terms of the generators A~ (A), S~(A), C~(A) and T>~(\) it is 
instructive to use the combinations 

V-(X) = 2XV-(X)-iA~(X), C'{X) = (2A + i)<T(A) (2.12) 

borrowed from the algebraic Bethe ansatz. Then the suggestive form of the quantum deter- 
minant reads 

(A~Z/)(A) = ^-(A + i/2)P"(A-i/2) - £~(A + i/2) <T(A - i/2) . (2.13) 

Thus in case of the c-number representation K(X — i/2, — ) connected to the left reflection 
algebra U~{X) the relation 

(A-K)(X - i/2, -) = 2(A - i) det K(X, -) = -2(A - i) — - -j^t^- — (2.14) 

[a y 

holds. Note that this connection is only valid for the shifted argument A — i/2 because the 
arising expressions in (12.111) are no longer of difference form. As the quantum determinant 
respects co-multiplication, applying it to the full representation (12.91) of the left reflection 
algebra with the monodromy T(A) yields 

(A-W)(A) = (d q T)(X - i/2) (A-K)(X - i/2, -) (d q T)(-X - i/2) . (2.15) 

The right reflection algebra can be treated in a similar way. We may leave with the 
suggestive form of the result 

(A+W)(A) =P + (A-i/2)^ + (A + i/2)-£ + (A-i/2)C+(A + i/2) . (2.16) 

Again, we used some suitable combinations reading 

A + (\) = -2\A + (\) - lD + (A) , C + (A) = (-2A + i)C + (A) (2.17) 

where in case of the c-number representation K(X+i/2, +) to the algebra U + (X) the quantum 
determinant takes the form 

(A+K)(X + i/2, +) = -2(A + i) det K(X, +) = 2(A + i) (A ~ °f " +) . (2.18) 

(a+) 
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2.2 Separation of Variables 

For diagonal boundary matrices the spectrum of the model has be obtained with the algebraic 
Bethe ansatz by action of the operators B~(\) on the completely polarized pseudo vacuum |0) 
|24j . In the case of generic boundary conditions, |0) is not an eigenstate of the transfer matrix 
and this approach is not possible. Instead one can follow Sklyanin's functional approach [25] 
based on the operator valued zeros Xj of i3~(A), j = 1, . . . , L. Studying representations of the 
reflection algebra on a space of symmetric functions of the eigenvalues Xj of these operators 
one obtains the TQ-equation for the eigenvalues A of the transfer matrix [9], 

A( Xj ) Q( Xj ) = i^A+(x j )Q(x j + i) + i^A-(x j )Q(x j - i) . (2.19) 

The coefficients A ± (A) factorize the quantum determinant to the transfer matrix (12 .7p ac- 
cording to (A+W)(A) (A~U)(X) = -A + (A-i/2)A~(A+i/2) . Generalizing the transfer matrix 
to include inhomogeneous shifts parametrized by L lattice parameters Sj the eigenvalues of 
Xj are found to be x^ = Sj ± i/2 and the functions A 1 * 1 are explicitly given as 



3 

\ + a - - L 



A- (A) — (A + i/2) (A + a +-./2)(A + a --,/2) n _ + + ^ + 

mini J- J. 



1=1 
L 



A+(A) = (A - i/2) (A - a+ + l/2 > (A - Q - + i/2) - s, - i/2)(A + S , - i/2) . 

rvT/-\/— J. J. 



(2.20) 



i=i 



Note that only the diagonal parameters a* 1 of the boundary matrices enter in these equations. 
To obtain the spectrum for non-diagonal boundary fields corresponding to values of the 
parameters ^ and 9 ± they have to be complemented with information on the asymptotic 
behaviour of A (A) at large |A| 3> 1: from the construction of the transfer matrix one easily 
obtains [9] 

a ( a) ch^wr+ff-'n . (2 . 21) 

a + a~ ch p+ ch f3~ 

Hence the parameter <fi is sufficient to characterize the influence of the non-diagonal boundary 
fields. As a change in the sign of one a-parameter can be absorbed into the corresponding 
j3 — > —(3 and 9 — > 9 ± in with the mapping ch <ft — > — ch <p the complete parameter range of 
(II. ip is governed by ch0 > along with the cases a + /i, a~/i > and a~ /i < < a + /i. As 
a simultaneous change a ± — > — a ± formally reverses all spatial directions we do not need to 
consider the range a + /i, a~ ji < 0. 

The xf are singular points of the difference equation ( 12. 19j) . i.e. simple roots of the 
coefficients A ± (x^ z ) = 0. Therefore, the Q-functions can be eliminated from (I2.19p in favour 
of a functional equation for A(x) valid on the discrete set g = {xf}- It has been shown [9] 
that this equation yields the complete spectrum of the transfer matrix for small lattices. In 
the homogeneous limit Sj — >• s it becomes 

a+a - a+a - 4s 2 + 1 
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Here he subscript g emphasizes that (I2.22p has been derived for the eigenvalue A (A) with 
arguments A taken from the set g and therefore holds only up to terms which vanish as 
(A — s) 2L+2 . Still, neglecting this fact and taking s to be a continuous variable one can solve 
for dlnAg(s) by Fourier transformation [9]. A particular state can be selected by imposing 
constraints on the analytical properties of A s . For example the ground state eigenvalue has 
no zeros in the strip | lmz\ < 1/2. Evaluating the result at the point i/2 (c.f. (14.24)) ) this 
leads to the correct bulk and boundary contribution E g to the ground state energy (ip(x) is 
the digamma function) 

E g = idlnA s (i/2) = + V(|a + |/2) - ^((|a + | + l)/2) + l/\a + \ 

+ V(|a-|/2) - V((|a-| + l)/2) + l/|a~| (2.23) 
+ 7r-21n2-l + (2-41n2)L . 

Finite size corrections of order 0(1/ L) which would capture correlations between the two 
ends of the chain are beyond this approach. 

2.3 Fusion Procedure 

The so-called fusion procedure grants the possibility to easily obtain i?-matrices and bound- 
ary matrices of higher dimensions obeying a Yang-Baxter or reflection equation respectively. 
In case of the i?-matrix this procedure is applicable to the auxiliary, the quantum space, 
and even both. Furthermore the associated transfer matrices are not independent from each 
other but satisfy functional relations called fusion hierarchies. 

An R- matrix of dimension k/2 in auxiliary space and a spin- 1/2 representation in quantum 
space is given by e.g. [TjJ] reading 

R {1 ... k)i {X) = P+. k R u {X)R M {X + i) • ■ ■ Rm{X + (k — l)i)P£. k . (2.24) 

Here, the projector P + is defined as 

P t-n = ^ £ Pa (2-25) 

where the sum runs over all permutations a G & n of the symmetric group and P a is the 
permutation operator reordering the positions in the space (C 2 ) ™ according to a. The fused 
boundary matrices are defined in [31] following Mezincescu and Nepomechie [T7| who carried 
out the first fusion step 

K+ 2) {\) = P+K+(X + i/2)i? 12 (-2A - 2i)K+(X - i/2)P+ (2.26) 

similar to (12.24)) explicitly. Utilizing the co-multiplication property one also finds 

W ( - 2) (A) = P+Ur(X)R 12 {2\)U;(X + i)P+ (2.27) 

yielding a transfer matrix with a spin-1 auxiliary space t^ 12 )(A) ee tr(i 2 ) Kt 2 > ( A) U7 l2 x ( A) . 
With these definitions Mezincescu and Nepomechie [17] showed the fusion formula for the 
transfer matrix of the open boundary model 

r (12) (A-i/2) = -(2A-i)(2A + i)r(A-i/2)r(A + i/2)-i(A+K)(A + i/2,+)(A-W)(A) (2.28) 



s 



depending on the quantum determinants of both reflection algebras and the original transfer 
matrix r(A) from ( 12. 71) . This result can be rewritten by absorbing the scalar pref actor in 
front of the r's into the definition of the fused transfer matrix giving 

t 2 (A-i) = ti(A-i)ti(A)-$(A) , r(A + i/2) = t^t!(A) (2.29) 

or a 

with the scalar function 5(A) on the RHS reading 

( A 2 + \)2L+1 

5(X)= { -^^-(X~a + )(X + a + )(X-a-)(X + a-) . (2.30) 



Extending this procedure to higher dimensional auxiliary spaces [20||3T] the arising trans- 
fer matrices t k for integer k are finally related to each other through the fusion hierarchy 

t fc (A-i(fc-l)) =t fc _a(A-i(A;-l))ti(A)-<y(A)t fc _2(A-i(A:-l)) , A; = 2,3,... (2.31) 

with S(X) given in (I2.30p and t (X) = 1. Note that we do not have to distinguish between 
the transfer matrices tk{X) and their eigenvalues because we are dealing with commuting 
quantities sharing a common system of eigenfunctions. Thus in the following we will use 
the notation 4 also for the eigenvalues. Again, the fusion hierarchy needs to be completed 
with the asymptotic behaviour 4(A) ~ akX^ 2Ij+2 ^ k of the eigenvalues of the fused transfer 
matrices: solving a recursion relation following from (I2.3ip and the asymptotic (I2.2ip . (12. 29ft 
of the eigenvalue A(A) of r(A) as a starting value one obtains 

M-a^W , nt = JL !fei^ . (2 . 32) 
2.4 Equivalence of TQ-Equations 

The fusion hierarchy (12.311) can formally be solved for the shifted eigenvalue ti(A— i/2) related 
to spin- 1/2 reading 

_ = t t (X ± i/2 - i(* + 1) + i) _ tjMX + i/2 ^ ■(* - 1) - p.33) 

n ' ' t fc _x(A + i/2 — ifc) v 1 ! tfc_i(A + i/2 — iA:) v ' 

matching the general form of a TQ-equation [29]. Indeed, renormalization of the fused 
transfer matrices 4(A) according to 

fe-i fc-i 
t k (X) = a k [H(X + i£) 2L+1 HiX + ii + iffl-^niX) (2.34) 
i=\ i=i 

and assuming the limit limfc^oo Tk(X — ik — i/2) = g(A) to exist yields a difference equation 

A(A) g(A) = ( -^$r (A - i/2) 2L+1 g(A + i) 

(—-\\L-<j> r , (2.35) 

+ ( A + i/2) 2i+1 [ J] (A - i/2 - « CT )(A - i/2 + a°)\ q(X - i) 

cr=± 
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which fixes the eigenvalue A(A) of t(A). With the transformation 

q(\) = (Ti) iA e iA ^g(A) [r(-i« + + 1/2 + iA)r(-i«- + 1/2 + iA)] _1 (2.36) 

we can absorb the exponential dependence on into Q(X) recovering the TQ-equation already 
obtained from the functional Bethe ansatz [9] reading 

±A(A) Q(A) = ( ~ 1 / (A - i/2) 2L+1 (A - a + + i/2) (A - a~ + i/2) Q(A + i) 



(2.37) 



2Aa+a:~ 

In addition to the explicit appearance of a 1 * 1 the boundary conditions enter this equation 
through the large-A behaviour (I2.2ip of the eigenvalue A(A) and r(A) respectively. 

For boundary parameters giving ch0 = ±1 equation (I2.37P can be solved by an even 
polynomial Q(X) = Ylj=i{^ ~~ A/)(^ + A?) where the Xj are roots determined from the Bethe 
equations 

/ Xj + i/2 \ 2L+1 X, +a + - i/2 Xj +oT- i/2 = X, -X k + i X, + X k + i 

\Xj - i/2 J Xj - a+ + i/2 \ - or + i/2 AJ- A 3 - - A fe - i Xj + X k - i ' [ ' J 

This case corresponds to boundary fields which can be dealt with by means of the algebraic 
Bethe ansatz. On the other hand, for ch0 ^ ±1 no simple ansatz for Q(X) is known. 
However, for large |A| it has to grow exponentially Q ~ exp(— i0A) according to (I2.36p . 



3 Truncated Bethe Equations 

After explicitly carrying out the large-fc limit of the auxiliary space dimension for the transfer 
matrix t k we are able to derive truncated Bethe equations in the non-diagonal boundary 
parameter range of <p G R. Studying the root distributions of t k for small system sizes yields 
information about the general root structure and it is possible to approximately treat the 
system with a finite number of zeros. With this knowledge on the analytical properties 
of the g-functions (12.351) leads to truncated Bethe equations. As a remark in the case of 
the XXZ model with diagonal boundaries and a special choice of boundary parameters the 
fusion hierarchy truncates exactly and the problem was already solved in [30] and [19] for 
non-diagonal boundaries respectively. 



3.1 Iteration of the Fusion Equations 

For systems with a few lattice sites (L < 14) the polynomial eigenvalues of the transfer 
matrix can be calculated explicitly either by exact formulas or numerical determination of 
the coefficients. Then the fusion hierarchy (I2.3ip can be used to compute the corresponding 
Tjfc(A) up to a specific fusion level k. 

The roots of these functions r k (X) are founcfj] to be in the strip —k — 1/2 < ImA < 1/2 
and are symmetric with respect to the line ImA — —k/2. As a consequence the zeros of 



1 This may not be true for all states. In cases where the 'Bethe ansatz-part' of the roots contains one or 
more string solutions the roots of q(X) may extend into the half plane Re A < 0. 
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Figure 1: Typical root distributions for the A-state of r k (X — ik — i/2) for an even system 
size L and Ima",Ima + > 1/2 (a) in the Bethe ansatz solvable case of ch0 = 1 and (b) in 
the non-diagonal case of ch0 > 1. Both distributions are symmetric to A = i(k + l)/2 . 

Tfc(A — ik — i/2) in ( 12.341) . which are relevant in the limit k — > oo for the q- function, emerge in 
the region < Im A < (k + l)/2 as denoted in Figured] for the Bethe ansatz solvable case (a) 
and for non-diagonal boundaries (b). Here we will concentrate on two distinguished states 
which we label by A and B: in the diagonal limit ch = ±1 amenable to the algebraic Bethe 
ansatz the first state A turns into the singlet ground state of the antiferromagnetic chain 
whereas the -B-state describes in this limit the fully magnetized state. Below we will study 
these states for boundary parameters Ima^ > 1/2 which excludes boundary bound states. 

For the Bethe ansatz solvable case in Figure [1^ we observe in the A-state a distribution 
of zeros on the real axis, which are the known 2 x L/2 Bethe roots for this sector. In addition 
the g-function has roots on the imaginary axis which form a half-infinite lattice of spacing i 
starting at the points A = a ± + i/2. In the limit k — > oo this lattice becomes exact which 
allows to rewrite ( I2.35P as (I2.37P with polynomial Q(X) as discussed above. The extra roots 
appearing on the symmetry line ImA = (k + l)/2 for finite k can be neglected in the limit 
k — > oo. 

In cases with no Bethe ansatz the roots on the symmetry line move to branches in the 
complex plane which strongly depend on the asymptotic ch0. A typical root configuration 
for the v4-state is shown in Figure [TJo: we find that the number of roots forming these 
branches is fixed and their positions in the complex plane show only a slow variation with 
respect to the system size. The additional roots near the real axis evolve into the 2 x L/2 
solutions of the Bethe equations (12.38}) parametrizing the antiferromagnetic ground state in 
the limit ch.<p — > 1. The zeros on the imaginary axis start from the points or* 1 + i/2 
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with exponential accuracy (in L) and are approximately spaced by i as in the Bethe ansatz 
solvable case. However, at some imaginary point depending on the boundary parameters 
there is a crossover of the sequence v k into a lattice of roots spaced by i/2. This lattice 
is approached with exponential accuracy which suggests to replace the sequence v k beyond 
some (half) integer position \M by a compensating factor accounting for the asymptotics. 
This allows to treat the finitely many remaining roots with imaginary part < M as in the 
treatise of the Bethe ansatz solvable case from above. 

For the 5-state a strong size dependence of the branches of extra roots is found which 
requires a more careful finite size analysis (see below). 



3.2 Truncation 

The analyticity of the transfer matrix eigenvalues (I2.35H implies that the roots q(Xj] 
have to satisfy an infinite hierarchy of Bethe equations 



c 



for all j £ N. In addition the restriction 

e"^(i/2 - a+)(i/2 + a + )(i/2 - a")(i/2 + a") = (3.2) 

emerges from the residue at A = to vanish. Our findings on the asymptotic positions of 
the roots with ImAj ^> 1 (c.f. Figure [Tb), i.e. by consecutive integers and half integers on 
the imaginary axis suggest us to explicitly deal only with the finite number N of zeros with 
ImAj < M for some sufficiently large M depending on the boundary parameters a ± and 0. 
Within this approach the roots with ImAj > M cancel (almost) perfectly in ( I2.35p . i.e. we 
can replace 

/ \ | >\ N N-2 

^ K 1 k=i e=i /„ „\ 

N-2 N 



q(X - i) 



^n(A-A,-i)n^ 



g(A) /(A-i)iV L\ x - X * ' 

The asymptotic behaviour (12.211) of the eigenvalue A(A) in the representation (12.351) requires 
/(A) = 1. As a result we obtain 

A(A) (A - i/2) 2L+1 (A - A. + i)(A - Ajy-i + i) ff 

era 2A - LJ - A — Afc 

w ( 3 - 4 ) 
(-l)^e-* 2i+ i n. =± (A-i/2-^)(A-i/2 + q-) A - A fc - i 

a+a-2A l+1/j (A — Ajv)(A - Xn-x) ^ A - A fc 
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from which we can derive truncated Bethe equations by the requirement of vanishing residues. 
The N remaining roots of q(X) have to satisfy the system 



-2<f> 



Xj - i/2 

N— 2 

(Xj — Aiv-i)(Aj — Aat_i + i)(Aj — Ajv)(Aj — A^ + i) t-t Aj — A& + i 



(3-5) 



( A . _ i/2 _ a +)( Aj - i/2 + a+) (Xj - i/2 - a~)(Xj - i/2 + a") A, - A, 

Mi 

for j = 1, . . . , N — 2. For j = N — 1, N the equations are ill-defined since the RHS vanishes 
identically. Note, however, that the analyticity of A(A) requires the residues of (j3.4p at 
A = A^v and A = Ajy_i to vanish. This condition can be met by setting Aat„i = iM + i/2 and 
A A? = iM + i for the (half) integer M chosen for the truncation above as the zeros Ajv-2 ~ i-^ 
and Aat„ 3 ~ iM — i/2 are located close but to not precisely on the asymptotic positions. As 
a test from the singularity of the eigenvalue A at A = the roots have to obey the restriction 



-2<j> _ Ajv(Aat — i)Ajv-i(Ajv-i — i) -rr A& 

" (i/2 - a+)(i/2 + «+)(i/2 - a-)(i/2 + a") Al A^ + i 



c 



• (3-6) 



k=l 



Equations (13. 5p can be solved numerically by Newton's algorithm. For the A-state this 
is possible due to the fact that the number and position of the roots forming the branch 
(c.f. Figure db) vary only slowly with respect to the system size for fixed boundaries. For 
the 2 x L/2 zeros near the real axis one can use the root distribution from the Bethe ansatz 
solvable case f l2.38D as starting values. Doing so one obtains solutions for up to several 
thousand lattice sites which can be used to examine e.g. the finite size effects (c.f. Figure H] 
in Section H]). 

In the .B-state the shape and position of the branches as well as the number of roots 
forming them strongly varies with the system size which prevents the derivation of truncated 
Bethe equations used above. However, a careful analysis of the arising root distributions 
leads to characteristic patterns as depicted in Figure [2] where we show typical results for 
a-parameters of equal sign. 



3.3 Scope of Application 

Clearly the method can also be applied to boundary bound states < Ima ± < 1/2 with 
roots sticking to the points i/2 — now being part of the parametrization (13. 4p of the 
eigenvalue. For a plain overview of the terminology and parameter range see e.g. [22] . 

Note that the derivation of the equations (13. 5 j) does not depend on the signs of the bound- 
ary parameters a^: for boundary fields with opposite signs, sign(Ima + ) = — sign(Ima:~), 
the root configuration corresponding to the A-state differs from the case above by the absence 
of the additional branches in the complex plane. For the 5-state, branches as in Figure [2] 
are still present. Therefore, in both cases an analysis of the truncated Bethe equations is 
possible along the lines described before. 
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Figure 2: Root distributions for the S-state at a ± = lli/3 and ch0 = 11. Zeros belonging 
to fixed lattice sizes L = 1 to L = 9 are each connected by a dashed line as a guidance for 
the eyes. The roots on the imaginary axis forming an asymptotic lattice of half integers are 
not displayed. 

4 Y-System and Non-Linear Integral Equations 

In order to directly capture the corrections to the bulk and boundary part ( I2.23|) of the 
energy eigenvalue of the open spin chain one can use the methods of complex calculus to 
utilize the equivalent representation 



t k (X + i/2)t k (X - i/2) = t k ^(X + i/2)t k+1 (X - i/2) + J] 5(X - i/2 + i£) (4.1) 

t=\ 

of the fusion hierarchy (I2.3ip . Following the standard scheme [TBj one introduces the combi- 
nation 

*M = %^ (4.2) 

being part of an infinite series {y k } related by functional relations. Once y<i is known the 
eigenvalue of the underlying integrable model can be calculated from the lowest level 

A(A + i/2) A (A - i/2) = (i + y2 ( A) ) (4.3) 

of the fusion hierarchy (14. ip e.g. by Fourier techniques. Note the invariance of the functional 
equation H4.3[) with respect to A — > —A or a ± — > ~a ± . 

4.1 Y-System and Fourier Transformation 

The aforementioned infinite system of functional relations is denoted as F-system and can be 
derived for products of the transfer matrix eigenvalues. It is sometimes called the universal 
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form of the TBA-equations [16] and follows directly from the fusion hierarchy (14. ip . By 
defining 

t fc _ 2 ( A-i(fc-2)/2)t fc (A- ifc/2) 

and explicitly calculating (1 + 2/fc+i) (1 + Vk-i) by making use of the fusion relation (14.11) one 
finds for integers k 



VkW = ^fc-i — ( 4 - 4 ) 



y k (X - i/2) y k (X + i/2) = (l + y*_i(A)) (l + y fc+ i(A)) (4.5) 
with jji = 0. Note that a simultaneous scaling of t k (X) and 5(A) from (12.311) 

tl (A)->^, (5(A) -> ; , ^(A-Klfc-l)) 7 (4.6) 

pM ' p(A)p(A-i) " ntoP(^-i<) 

with any function p(A) leaves the F-system (14.51) invariant. In view of this fact we choose 
p(X) = (2A + i) _1 which allows to consider polynomial 5(A) and t k (X) rather than rational 
functions. Thus the zeros and poles of the ^-functions can be identified with the zeros of 
ifc(A) and S(X), respectively. Taking the logarithmic derivative of the F-system (14 .5p and 
Fourier transforming 



oc 



f(k)= jdxe- ikx f(x) , f(x)= J^-e ikx f(k) (4.7) 

— oo — oo 

the imaginary shifts occurring on the LHS can be removed leading to an infinite set of non- 
linear integral equations (NLIE) for the y-functions evaluated on the real line |27j . 

\ny 2 (x) = d 2 (x) + s* ln(l + y 3 ) , ^ 
In y k (x) = d k (x) + s * ln(l + y fc _i) + s * ln(l + y k+1 ) , k > 2 . 

Here (s * f) denotes a convolution of / with the kernel 

six) = —4 — r , s(k) = - } ; r . (4.9) 
v ; 2ch(vre) ' v ; 2ch(A;/2) v ; 

As a consequence of the zeros and poles of the functions y k their logarithmic derivatives are 
not analytic which produces the additional contributions dj(x) due to the residue theorem. 
The determination of these model-dependent driving terms is the challenging component of 
the problem. To reduce these terms we consider the toy equation 

y(x + i/2)y(x-i/2) = F(x) (4.10) 

for y(x) where F(x) is given explicitly with constant asymptotic and the auxiliary condition 
y(±x) = for some \ £ { z £ C O < Imz < 1/2}. Taking the logarithmic derivative one is 
led to 

oo 

Jdxe-' lkx dlny(x±i/2) = e Tk/2 dhuj(k) t 2nie Tikx . (4.11) 
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Switching to Fourier space we can solve (I4.10p for d\ny(k) using the integral representation 
of the digamma function yielding 



dfc sh(i/fc/2) ^ = _ J_ 
2tt ch(fc/2) 2m 



d T In 



ch(7rx) — sin(7ri//2) 
ch(7rx) + sin(7r^/2) 



(4.12) 



Identifying u — 1 — 2%/i, inverse Fourier transforming and integration with respect to the 
variable x yields 



lny(x) = C + In 



ch(7rx) — cos(i7rx) 
ch(7ra;) + cos(i7rx) 



+ (s * \nF)(x) 



(4.13) 



where the integration constant C can be recovered from the asymptotic condition of the 
initial equation (I4.10p . 

Turning to the actual zeros of the ^-functions we find that each yk>2 has a double zero 
at A = whereas 2/2 exhibits a (2L + 2)-fold zero at A = 0. To proceed we have to choose a 
particular state for which the transfer matrix t\ and yk do not contain further zeros in the 
strip I Imz\ < 1/2 usually called hole-type solutions. Such a state exists for even lattice sites 
and opposite signal of the boundary parameters a 1 * 1 differing in their absolute values only by 
some finite amount (see below). This state is restricted to G M, i.e. ch0 > 1, and in the 
diagonal limit of — > it becomes the state with lowest energy in the sector of vanishing 
magnetization labeled A-state before. The corresponding driving terms (I4.13P read for each 
double zero at x = with x — ► 



In 



ch(7rx) — cos(i7rx) 
ch(7ra;) + cos(i7rx) 



lntanh 2 \nx/2\= In 6? (a) 



(4.14) 



The poles of the ^/-functions can be treated in a similar way. Again, using (14. 6[) to scale the 
transfer matrix properly only the polynomial denominators, e.g. 5(A + i/2)<5(A — i/2) for 7/3 are 
responsible for the pole structure with respect to the boundary fields a ± . The multiple poles 
at A = are compensated by the transfer matrices. Already from the form of the arguments 
of each S given in ( 12 .30 j) it is clear that the positions of poles vary by i/2 for successive k and 
two poles of the same k differ by i. The general scheme is depicted in Figure [3j Thus using 
the parametrization a = i(n/2 + a') with n e No fixed and < a' < 1/2 the function |//t(A) 
shows poles at A = ±ia' for indices k e {n + 2N} and poles at A = ±(i/2 — ia') for indices 
k G {n + 1 + 2N}. The corresponding driving terms arise from (14. 12ft by inverting the sign 
on the RHS accounting for poles rather than zeros yielding negative logarithms of 



x 



K(x) 



ch(7ra) — cos(7ra') 
ch(7rx) + cos(7ra; / ) 
ch(7rx) — sin(7TQ; / ) 
ch(7ra;) + sin(7ro; / ) 



for k G {n + 2N} 
for k G {n + 1 + 2N} 



(4.15) 
(4.16) 



2 Note that for a ± -parameters of equal sign and small k the branches in Figure [T}) appear as unwanted 
hole- type solutions in the considered strip. The same holds in the parameter range <f> G ilR for any combination 
of the a-parameters. 
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Figure 3: Exemplary position of poles (o) of yt in the complex plane arising from one 
boundary field a for a) k — 2, b) k — n + 1, c) k — n + 2 and d) A; = n + 3. With 
these examples all different cases of singularities inside the strip | ImA| < 1/2 are already 
considered. 



Note that only the modulus of the parameter a enters and the terms h c and h s alternate 
for successive k. In summary a single zero at A = gives a driving term of lntanh \irx/2\ 
and due to logarithmic derivative the multiplicity is only reflected in the integer prefactor. 
Rearranging all driving terms in matrix form 



b n k (x) 





1 


2 


3 


4 


" M 





tanh 2 \tcx/2\ 


h c (x) 


h B (x) 


h c (x) ■ 




1 


1 


1 


h c (x) 


h s (x) ■ 




2 


1 


1 


1 


h c (x) ■ 




n 










/ 



(4.17) 



V 

allows us to write the non-linear integral equations in a compact form. For boundary pa- 
rameters = ±i(n ± /2 + a' ± ) with n + = n~ = n e Nq we find for an asymptotic condition 
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ch > 1 of the transfer matrix the infinite set of equations 



lny 2 {x) 



(L + l)\nb° 1 (x)-\nbf(x)-\nbf(x) + s*\n(l + y 3 ) 

lnbl(x)-lnhf(x)-lnb%~(x) 

+ s* ln(l + + s * m(l + y k+1 ) , k > 2 



(4.18) 



for an even L. With the asymptotic condition b^(x) — > 1 for x — > ±oo no further integration 
constants have to be considered as the constant asymptotic 



already satisfies the hierarchy (14. 1 8j) in the limit of x — > ±oo. 

To solve an infinite system of equations such as (I4.18P numerically a controlled scheme 
for its truncation is needed. In the standard TBA approach [28j, e.g. for periodic boundary 
conditions, such an approximation can be based on the fact that (i) only the first few of 
the NLIE contain a non-zero driving term and (ii) the asymptotic solution (14. 19ft solves the 
NLIE without driving terms. Replacing yk{x) by this constant asymptotic for some k chosen 
sufficiently large one is left with a finite set of NLIE which can be solved numerically very 
efficiently 

In the case of open boundary conditions the NLIE contain driving terms for all levels k 
and therefore the asymptotics (14.191) determine only the large-x behaviour. Nevertheless, our 
numerical results using the constant asymptotic of as an approximative limiting function 
seems to assure convergence of the system reasonably well. In the special cases of diagonal 
boundaries <fi = we find that it is better to scale the constant asymptotic by the necessary 
driving terms in order to gain higher accuracy or use less equations without losing accuracy. 

For asymptotic parameters <fi ^ the constant solution (14. 19j) grows exponentially with 
k and numerical limitations are quickly reached restricting the method to ch</> = 0(1). The 
number of non-linear integral equations, necessary for decent accuracy, produces function 
values for which more sophisticated numerical treatment is needed. In principle, this simple 
approach which just uses the asymptotic as a limit function can be improved by considering 
solutions of (14.181) in the limit of max{n + ,n~} <C k. It turns out, however, that this more 
sophisticated treatment of the asymptotics shown in Appendix |A] does not allow to reduce 
the number of NLIE substantially and that the restriction to ch0 = 0(1) prevails. 

4.2 Energy Eigenvalue 

Once the full hierarchy for yk{x) is solved the function y 2 (x) (and its analytical continuation) 
can be used to determine the polynomial eigenvalue A (A) of the transfer matrix. Identifying 
the quantum determinant 5(\) = (a + a~) 2 A g (\ + i/2)A g (A — i/2) with (12.221) we see from the 
lowest level 



of the fusion hierarchy that the corrections to bulk and boundary contributions [9| are con- 
tained in y 2 . The logarithmic derivative of the eigenvalue reads after simple manipulations 



,oo 




(4.19) 



sh 2 



A(A + i/2) A(A - i/2) = A 9 (A + i/2)A s (A - i/2) (l + y 2 (A)) 



(4.20) 
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in Fourier space 

oo 

/r\ Ir pikx 
^ 2ch(k/2) dHl + y2)ik) (421) 

— oo 
oo 

= 8toA,(s)-,/f y;-'>» Ml + . (4.22) 

J I ch (tt(x - y)) 

— oo 

The expression can be used for asymptotics > 1 (i.e. ch0 > 1) and opposite signs of the 
boundary fields = ±i(n ± /2 + a ±r ) with < a ±! < 1/2 and n ± E N fixedly However, a 
change of the sign of one a-parameter transforms — > <fi ± m but leaves the asymptotics y%° 
from (14.191) unchanged. This again reflects the invariance of the functional equation (14 . 2 [) 
with respect to A — > —A and or* 1 — > —a ± . As a remark the polynomial eigenvalue A(A) can 
be evaluated explicitly in the limit of — > oo reading 

A(A) = ti^(A 2 + l/4) i+1 . (4.23) 

Note that the prefactor is recovered from the asymptotic behaviour (I2.21f) and does not 
contain the combined limit of <f), \a^\ — > oo but ch<f)/(a + a~) fixed. Continuing (I4.2ip and 
(I4.22p respectively to x = i/2 leaves us with the energy eigenvalue of the hamiltonian (II. ip 

oo 

E = i«91nA(i/2) = i«91nA 9 (i/2) - rr J ln(l + y 2 )(y) • (4.24) 

— oo 

The bulk and boundary part E g = i91nA s (i/2) was already specified in (12.231) and re- 
mains valid for any boundary conditions and all < |a | < oo as the poles emerging from 
< (o^l < 1/2 in both A(A) and A S (A) cancel out. Corrections due to correlations between 
the boundaries in the finite system are covered by the second term yielding the energy de- 
pendence in Figure HI Clearly in the limit — > oo the expression (14.231) diverges when the 
energy is calculated and reveals all hole- type solutions to accumulate at A = ±i/2. 



5 Summary and Outlook 

In this paper we have presented two approaches for the analysis of the functional equation 
describing the spectrum of the XXX spin chain with non-diagonal open boundary conditions 
(II. ip . Based on the integrable structures underlying this model the TQ-equation arise both 
from Sklyanin's separation of variables and from the fusion procedure for transfer matrices. 
Usually, the solution of the TQ-equation for lattice models with compact realization of the 
symmetry can be re-expressed in terms of algebraic Bethe equations which are obtained in 

3 In our derivation of the non-linear integral equations (14. 18)) we have restricted the boundary fields by 
choosing n + = nT . In numerical solutions of the equations we find, however, that the error due to the 
violation of this constraint is of similar order as the one arising from the truncation of the infinite hierarchy, 
see Figure HI 
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Figure 4: Influence of non-diagonal boundaries on the finite size effects to the energy eigen- 
value for the A-state with parameters a + = 5i/3 and a~ = — 7i/5 of different sign and lattice 
sizes L = 2,4, 10,20,50. For large asymptotics ch(/> the calculations are performed with the 
truncated Bethe equations (o) matching for intermediate values of ch0 the results from the 
F-system (+). The latter one deals with low asymptotics down to ch</> = 1. Solid lines are 
data from exact diagonalization for small system sizes L — 2,4, 10. 

a straightforward way starting with a polynomial ansatz for the Q-functions. For generic 
boundary conditions such an ansatz does not satisfy the requirements on the asymptotic 
behaviour of the eigenvalues in the present case (see [HE]). Starting from a finite size study 
of the analytic behaviour of the transfer matrix eigenvalues we have been able to express 
the functional equation in a way which opens the possibility of a numerical treatment for 
systems sizes beyond what is accessible to exact diagonalization. Both the truncated Bethe 
equations derived in Section [3] and the non-linear integral equations presented in Section [5] 
can be applied to compute the spectrum for non-diagonal boundary fields parametrized by 
chef) > 1. They are particularly useful to compute the energy of the state which evolves 
into the antiferromagnetic ground state of the chain for diagonal boundary conditions (the 
A-state). 

The truncated Bethe equations are found to work especially well for ch0 3> 1 and can 
deal with arbitrary boundary parameters a^. Going to large system sizes the numerical 
analysis of the equations is limited in principle by the necessity to find appropriate starting 
values for the Bethe roots. Still, selected energies can be computed for systems of several 
thousand sites. In the second method of non-linear integral equations the system size enters 
merely as a parameter. At the same time, however, boundary parameters have to be 
of different sign to meet certain analyticity requirements. In addition, the presence of zeros 
together with the exponential growth of the asymptotic value of the y-functions with the 
level k leads to numerical instabilities limiting the use of these equations to boundary fields 
with ch.(p = 0(1). Taken together the two methods are complementary allowing to cover 
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the entire range of boundary parameters with ch0 > 1. Our numerical data show these 
approaches agree for intermediate values of ch0. Further support for their validity comes 
from comparison with exact diagonalization for lattices with up to 14 sites. 

At the same time there remain several open problems which we shall address in the 
future: for the XXX chain (11.11) considered in this paper the description of the root dis- 
tribution for the state evolving into the fully polarized ferromagnetic state in the limit of 
diagonal boundary conditions (5-state, see Figure [2]) needs to be modified to extend the 
non-linear integral equations approach to this state. The limitation of the non-linear inte- 
gral equations for the A-state to boundary parameters satisfying ch0 > 1 appears to be a 
technical problem which could be resolved eventually by finding an exact truncation of the 
infinite hierarchy of equations similar to that for the s/(2) model [26j. Furthermore, the case 
< ch0 < 1 corresponding to Hermitian boundary terms in (11 .ip is not covered by our 
present approaches. An extension to this range of parameters may be possible following the 
treatment of excited states for the periodic XXZ chain to deal with hole-type solutions in 
the complex strip to be considered for the Y"-System [13J. Similarly, the present restriction 
to opposite signs of the a-parameters in the non-linear integral equation approach could be 
resolved. Finally, the methods introduced here should be extended to the XXZ chain with 
non-diagonal boundaries. The corresponding TQ-equations have already been derived from 
the fusion procedure [29] . Of particular interest in this model is the expansion of the en- 
ergy eigenvalues around the Bethe ansatz solvable case, ch0 = 1, which would allow for a 
non-perturbative study of current fluctuations in certain models for diffusion in one dimen- 
sion [SHE]- 
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A Asymptotic Truncation 



As the system of equations (I4.18P has alternating driving terms we assume two limiting 
functions gi(x) and g 2 (x) to truncate the system by 

^t:" r :: !,Ar) , ^ ^tr::: !h[r] ca.d 



bf (x)b--(x) 



W{x)bi;{x) 



for successive k and k! = k + 1 ^> max{n + ,n } supported by numerical observations. In- 
serting these into (14.181) we obtain two coupled non-linear integral equations 



In gi(x) = s * In 



bl92 



( 1 + i/£?_i 



1 + 



/ 



°k+l"k+l 



6? 91 



1 + 



6? Sl 



A 



(A.2) 
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In the limit of k — >■ oo the coupled system linearizes due to the exponentially fast growing 
asymptotics into 

In g 1 (x) = 2s * In g 2 + 2s* In b\ - 2s * In ftjj^ - 2s * In ftj^ g ^ 
In <72(^) = 2s * In + 2s * In 6° — 2s * hi 6^1 : — 2s * lnfe^,_ 1 

and is solvable in Fourier space after differentiating to apply the residue theorem. The explicit 
expression reads 



x ch(7ra;/2) ch(7rx/2) 
2ch(7re) 



df(x) df(x) 



2 



, 92{x) 



, . . ch(7ra;/2) ch(7ra;/2) 
2ch(7ra;) 



df +l {x) dl +l {x) 



2 



(A.4) 



leaving the asymptotic of yk{x) unchanged as #1/2 (a?) — > 1 for x — > ±00 with the shorthands 

™ / s I ch(7Tx) + cos^c/) for n + fc even . , , 

dl(x) = I ) ' ; ;x 7 . (A.5) 

I ch(7nr) + sin(7ra: ) for n + k odd 
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